WKB theory of epidemic fade-out in stochastic populations 
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Stochastic effects may cause 'fade-out' of an infectious disease in a population immediately after 
an epidemic outbreak. We evaluate the epidemic fade-out probability by a WKB method and find 
that the most probable path to extinction of the disease comes from an instanton-like orbit in the 
phase space of an underlying Hamiltonian flow. 
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An infectious disease can disappear from a population 
immediately after a major infection outbreak [l|,l2|- This 
phenomenon, called "epidemic fade-out" , occurs if the 
epidemic dynamics is oscillatory, and the number of in- 
fected individuals at the end of the first outbreak of the 
disease is relatively low so that fluctuations in the disease 
transmission can "switch ofP' the disease. Epidemic fade- 
out has been addressed by epidemiologists via stochastic 
simulations. One exception is Ref. ^] which arrived at 
important analytical results which we briefly review and 
use below. 

Epidemic fade-out is an example of a large fluctuation 
in a multivariate stochastic system far from equilibrium. 
There is no general theory of fluctuations in such systems, 
and finding the probability of a large fiuctuation is hard. 
Here we develop a theoretical framework for analysis of 
epidemic fade-out of the example of stochastic SI model: 
a Markov process involving Susceptible and Infected sub- 
populations fl|, 12) S|- We formulate the problem in a 
master equation setting. In contrast to endemic fade- 
out, which can be studied assuming quasi-stationarity 
Q , the epidemic fade-out occurs on a fast time scale (de- 
termined by the deterministic rate equations of the SIR 
model) , so no ready-to-use methods of solution are avail- 
able. To overcome this difficulty we derive a stationary 
equation for the mean residence time of the population 
in a certain state. Then we develop a WKB theory using 
the population size as a large parameter. In the WKB 
framework the problem reduces to that of finding a spe- 
cial zero-energy phase orbit of the underlying Hamilto- 
nian: the orbit which provides a global minimum to the 
action functional for boundary conditions, correspond- 
ing to epidemic fade-out. This orbit, which encodes the 
most probable path of the population towards the dis- 
ease extinction, turns out to be instanton-like We 
observe that the epidemic fade-out instanton exists only 
in the regime when the epidemic dynamics, as described 
by the deterministic rate equations, is oscillatory. Of spe- 
cial interest is the regime where the number of infected 
exhibits large oscillations prior to reaching the endemic 
state. By using a matched asymptotic expansion, we an- 
alytically calculate the action along the instanton which 
determines the epidemic fade-out probability. 



The probability Pn,m{t) to observe, at time t, n suscep- 
tible and m infected individuals is governed by the mas- 
ter equation with transition rates from Table 1. Solving 
the master equation analytically is hard. The often used 
van Kampen system-size expansion (vKSSE) ^ approxi- 
mates the master equation by a Fokker-Planck equation. 
In the context of epidemic fade-out in the SI model the 
vKSSE was employed in Ref. (iL The vKSSE is very 
useful for "typical" fluctuations but it often fails for 
large fluctuations 0], such as those causing extinction. 



Event 


Type of transition 


Rate 


Infection 
Renewal of susceptible 

Removal of infected 
Removal of susceptible 


I^I + l 
S ^ S+1 
I ^ I-l 
S^S-1 


{I3/N)SI 

ri 



TABLE I: Stochastic SI model 

Before dealing with the master equation, consider the 
deterministic rate equations for the SI model: 

S^nN~{/3/N)SI-fiS, i ={I3/N)SI ~TI . (1) 

For /3 > r there is an attracting flxed point S — {T/f3)N, 
I — /i(l/r — 1 /P)N which describes an endemic infection 
level, and an unstable (saddle) point S — N, I = which 
describes an infection-free population. At /i > 4 (/3 — 
r)(r//3)^ the attracting fixed point is a stable node. We 
are mostly interested in the opposite inequality, when the 
attracting fixed point is a stable focus, and the epidemic 
dynamics is oscillatory. Let a few infected be introduced 
into a susceptible population. For small /i the minimum 
number of infected at the end of the first outbreak of the 
disease is small, see the dashed line in Fig. [TJ As a result, 
stochasticity, missed by the rate equations, can "switch 
off' the disease before the endemic level is reached. To 
describe this process we use the master equation 

n' . m' 

= /i [iV(P„_i,„- Pn,ni) + (n-|-l)P„+i_,„- nPn^yyi] 

+ T[{m + l)P„,TO+i - niPn^yn] 

-f {I3/N)[{n + l)(m - l)P„+i,™_i - nmP„,„] .(2) 
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A natural initial condition is a product of Kronecker 
deltas: Pn.mit = 0) = 6n,NSm,ma- One boundary con- 
dition reflects the fact that m = is, for any n, an ab- 
sorbing state. Being interested in epidemic fade-out, we 
exclude from consideration all stochastic trajectories that 
do not reach the extinction boundary m — immediately 
after the first outbreak and leave the region of small m. 
This is achieved by introducing an artificial absorbing 
boundary that will be specified shortly. 
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FIG. 1: (color online). An epidemic outbreak in the SI model. 
Shown is the rescaled number of infected y = I/N versus 
rescaled time fit. Dashed line: prediction from the rate equa- 
tions ([5]). Solid line: the epidemic fade-out instanton. The 
rescaled parameters K = (3/ fj, = 30 and 5 = 1 — V / (3 = 0.5. 

The disease can only disappear from the population 
via transition from a state (n, 1) to the state (n, 0). Con- 

oo 

sider the mean residence time T„^m = J Pn,m(t) dt of the 



system in the state {71,111), where m > 0. The accu- 
mulated extinction probability Vn from the state (n, 1) 
is Vn = r?ra,i, and the total extinction probability is 
V — X^n^"- Integrating Eq. ([2]) over t from to 00 
and using the equality Pn.mit = 00) = and the initial 
condition, we obtain a stationary equation for Tn,m>o'- 

^ ^ ^^n,m; n' ,m' Pn' ,ni' ~t~ <^n,no ^m,7no ^ • (3) 
n',m'>0 

We assume throughout this work that N ^ 1. Here the 
stochasticity is weak (but very important), and Eq. ([3]) 
can be approximately solved by the WKB ansatz Tn^m = 
a{x, y) e-'^S(x,v) where a and S are smooth functions 
of the continuous variables x = n/N — 1 and y = m/N . 

In the leading WKB order one arrives at a stationary 
Hamilton- Jacobi equation H{x, y, dxS, dyS) = 0, where 

H{x,y,px,Py) = eP== - 1 + (1 + x) (e-P- - l) 
+Kil - S)y {e-Py - 1) + Ky{l + x) (eP""?'- - l) , (4) 

and we have introduced rescaled parameters 6 = 1 — r//3 
and K = P/ fj, and rescaled time by the rate constant 
/i The four-dimensional (4d) phase space, defined 
by the Hamiltonian yields an instructive visualiza- 
tion of the most probable path of the disease toward 
fade-out. As H is independent of time, H{x,y,Px,Py) = 



E = const. Furthermore, in view of stationarity of the 
Hamilton- Jacobi equation, we only need to deal with 
zero-energy orbits, E = Q. The simplest among them 
are fluctuationless orbits lying in the plane Px = Py = 0. 
These are described by the equations 

X = -X - K y{l + x) , y = -K{1 -S)y + Ky{l + x) , 

(5) 

which coincide with the (rescaled) rate equations ((1]). 
Disease fade-out demands a fluctuational orbit, for which 
the momenta Px and Py are nonzero. Before dealing with 
such orbits, consider the fixed points of the zero-energy 
Hamiltonian. There are exactly three such points, all of 
them 4d saddles Q. Two of them, B = [0, 0, 0, ln(l - 6)] 
and C — [0,0,0,0], describe infection-free steady states. 
Point C is fluctuationless: it corresponds to the saddle 
point of the rate equations. Point B is fluctuational, 
as its Py ^ 0. Finally, the fluctuationless fixed point 
A = [—6, (S/ K){l — S)~^ , 0, 0] corresponds to the endemic 
fixed point of the rate equations. 

Let one or few infected be introduced into an infection- 
free population. In the leading WKB order this initial 
condition can correspond to different phase-space points 
whose projections on the xy-plane are very close to the 
fiuctuationless fixed point C. Each of these phase-space 
points generates an orbit which exits the fixed point C 
along the manifold spanned by it two unstable eigenvec- 
tors. For epidemic fade-out to occur, such an orbit must 
reach the extinction hyperplane y = before crossing, 
say, the hyperplane y = —{x/K){l — 6)^^, —S<x<0 
(which is a 4d extension of the artificial absorbing bound- 
ary mentioned above). One can prove that, among all 
such orbits, the one providing the global minimum to 
the action S{x,0) (and therefore the global maximum 
to the fade-out probability Vn) ends in the fluctuational 
fixed point B. As a result, maxPn = Vn- Therefore, at 
N ^ 1, the epidemic fade-out problem reduces to that of 
finding a heteroclinic orbit going from C to B. We found 
numerically that such a heteroclinic orbit CB exists if 
and only if if > Kc = (1/4(5)(1 — 15)~^: when the endemic 
fixed point, predicted by the rate equations, is a focus. As 
K exceeds Kc, the heteroclinic orbit emerges via a global 
bifurcation. In fact, one finds multiple heteroclinic orbits 
aX K > Kc. They can be classified by whether their xy- 
projections exhibit a single loop, two loops, three loops, 
etc. A single- loop orbit, see Fig. [21 corresponds to dis- 
ease fade-out immediately after the first outbreak. A 
two-loop orbit corresponds to a fade-out immediately af- 
ter the second outbreak, etc. The connection between the 
rapid epidemic fade-out in a stochastic population and 
a zero-energy instanton-like orbit of an effective Hamil- 
tonian is a central result of our work. We stress that 
previous studies, which found instanton-like orbits in the 
context of population extinction, dealt with extinction 
from long-lived metastable states. 

How does the epidemic fade-out instanton look like at 
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FIG. 2: (color online), a: An epidemic outbreak on the xy- 
plane as predicted by the rate equations ([5} (dashed line) and 
the epidemic fade-out instanton (thick solid line) . Also shown 
is the endemic fade-out instanton 4] (thin solid line). 6: py 
(inset: px) vs. t for the epidemic fade-out instanton. The 
rescaled parameters are K — 30 and S — 1/2. 
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FIG. 3: (color onhne). Same as in Fig.[2]but for K = 1.7875 
and 5 = 1/3, so Kc — 1.6875. The endemic fade-out instanton 
U is shown by the dash-dotted line. Inset: a blowup near the 
endemic fixed point A. 

different parameters? For KS ^ 1 the fraction of in- 
fected versus time, y{t), first rapidly grows and becomes 
large and then falls down to a small value [see Fig. [U 
solid line] , closely following the prediction from the rate 
equations. Then y{t) strongly deviates from the deter- 
ministic path and rapidly goes to zero. The x-, Px- and 
Pj,-dynamics for the same values of K and 5 are depicted 
in Fig. [21 On can see that a rapid deviation from the 
deterministic path occurs around x = —5. Notably, \px\ 
remains much smaller than unity everywhere. \py\, how- 
ever, is steadily growing and, at non-small 5, becomes 
0(1) as the instanton approaches the fluctuational fixed 
point B. As a result, the vKSSE is invalid for most of 
the small-y region where disease extinction occurs. 

Near the bifurcation, {) < K — Kc <C A'c, our nu- 
merics reveals an intimate relation between the epidemic 
fade-out instanton and two other zero-energy heteroclinic 



orbits. The first is the deterministic orbit which lies in 
the ccy-plane and goes from C to A. The second is the en- 
demic fade-out instanton: a heteroclinic orbit which goes 
from A to B and describes stochastic endemic fade-out 
The xy-projection of the epidemic fade-out instan- 
ton initially closely follows the deterministic orbit CA, 
see Fig. [3l The momenta Px and Py are very small here. 
They slowly build up and become important only when 
the xy-projection of the instanton reaches a close vicinity 
of the endemic point A. Here the projection of the epi- 
demic fade-out instanton departs from the deterministic 
orbit (see the inset of Fig. [3|) and rapidly approaches the 
projection of the endemic fade-out instanton. 

To evaluate V ^ Vn in the leading WKB order, we 
need to calculate the accumulated action Sq along the in- 
stanton. In the rest of this communication we will focus 
on the important regime of K5 3> 1, when the fade-out 
probability can indeed be significant. It turns out that 
the presence of the small parameter {ICS)"^ enables one 
to find the instanton, and calculate iSo, analytically. An 
immediate simplification comes from the fact that the 
fluctuations of the number of susceptibles are negligible 
everywhere, so we can Taylor-expand the Hamiltonian 
(U) in Pa: 1 and truncate the expansion at first or- 
der. Another simplification employs the strong inequal- 
ity y <^ S which holds in the whole region where the 
fade-out instanton significantly deviates from the deter- 
ministic orbit. A complete calculation of the instanton 
will be reported elsewhere. Here we will analytically cal- 
culate So . As can be verified a posteriori, the main con- 
tribution to (So comes from a narrow region |a; -|- (5| ^ S, 
where the instanton rapidly departs from the determin- 
istic orbit. Furthermore, ^ 1 in this narrow region, 
so one can Taylor-expand Eq. (U in py and truncate the 
expansion at p^. Neglecting small terms, we can reduce 
the Hamiltonian (|4]) to 

H{x,y,Px,Py) ~5px + Kypy [x + 5 + {I - 5)py\ . (6) 

The reduced problem is integrable. There is no need in 
the full solution, however, if one only needs to evaluate 
Sq. The Hamilton's equation for x yields x{t) = 5{t— 1), 
where the arbitrary constant is fixed by choosing x{t — 
0) = —S. The Hamilton's equation for py reads 

Py = -Kpy [x + S+il- S)Py] . (7) 

Plugging here x — S{t — 1), we obtain an exactly soluble 
equation for Py(t). The boundedness of Py{t) fixes the 
integration constant, and we obtain 

Now let us calculate S along the instanton: S = PxX + 
Pyy = H + K{1 — S)ypy = J^, where we have used H = 
and denoted T = K{1 — S)ypy. Using the Hamilton's 
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equations, we observe that T{t) obeys the equation T = 
—K{x + 5) J- = —Kdt!F. Integration yields 

= K{1 - 5)v{t)pl{t) = Cexp(-A'(5iV2) , (9) 
where C — const. Therefore, S = C exp{—K6t^ /2), and 



Sdt ^ C 




(10) 



What is left is to find C. Importantly, the deterministic 
solution still holds in the region of —x — (5 -C 5 (or — i ^ 
1). For K6 ^ 1 the deterministic solution was found 
by van Herwaarden In the region of —x — (5 <C (5, 
Eqs. (3.25 a-d) of van Herwaarden can be simplified and 
rewritten, in our notation, as 



yit) 

Vm 



yra exp(X,5tV2) 



(11) 



ym{K, 5) 



KS 



l + x„i \ S 
X exp [K{x^ + 5)-{l + x-i) Qi{xm)] , (12) 



where Xm = Xm{6) is the negative root of the equation 
Xm = (1 - ^) ln(l + Xm), and Qi{xm) is given by 



s(s + 6) 



(1 + s)2 [s - (1 - ,5) ln(l + s)] 

ds. (13) 



(1 + Xrn) {s - Xm) 



(For (5^0 one obtains Xm{^) — —2(5 and Qi ~ —4(5.) In 
the region of {KS)^^^^ <C — .t — 5 <C (5 Eq. ([8]) becomes 

Py{t) = -(1 - 6)-^[6/{2ttK)]^/^ expi-K6ty2) . (14) 

Using Eqs. ([9]), (fTTj) and (|14p in their joint validity region 
(KS)-^/^ -^-x-S-^S, we obtain C = ymS/[2Tr{l-S)]. 
Putting everything together, we obtain the leading-order 
WKB result for the epidemic fade-out probability: V ~ 
exp(— A^5o), where Sq is given by Eq. pU]) and ?/,„ is given 
by Eqs. p2)) and p3|) . Note that So is exponentially small 
in KS ^ 1, so the WKB result holds only for very large 
N: NSq ^ 1. In Fig. 2] our results for Sq are compared 
with those obtained by a numerical integration of the 
Hamilton's equations. For large KS the agreement is very 
good. For NSq < 1 the epidemic fade-out probability is 
large. 

That truncation of H at Py yields, at KS ^ 1, an 
accurate leading-order result for So justifies the validity 
of the vKSSE for calculating Sq. Indeed, our leading- 
order result for V coincides with that obtained, by an 
entirely different method, by van Herwaarden [3| whose 
starting point was the vKSSE. We reiterate, however, 
that the vKSSE is invalid in most of the small-y region, 
while the full WKB Hamiltonian (jH) holds there. Only at 
(5 <C 1, when \py \ <C 1 on the whole instanton, the vKSSE 
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FIG. 4: Natural logarithm of the action So along the instanton 
versus KS at different S as predicted by our asymptotic theory 
at KS ^ 1 (lines), and found by a numerical integration of 
the full Hamilton's equations (circles) for S = 0.5). 



describes the instanton correctly. In this case one obtains 
So = (2,5V7re4i^)i/2(e/2)-^^ 

In summary, we have shown that rapid epidemic fade- 
out in stochastic populations is amenable to an accurate 
analysis via a WKB theory. We calculated the fade-out 
probability and established an unexpected connection be- 
tween the rapid epidemic fade-out and an instanton-like 
orbit of an underlying Hamiltonian. The fade-out instan- 
ton should be observable in stochastic simulations of, and 
actual data on, epidemics in small communities. 
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